Bright-Dark Soliton Complexes in Spinor Bose-Einstein Condensates 
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We present bright-dark vector solitons in quasi-one-dimensional spinor (F — 1) Bose-Einstein 
condensates. Using a multiscale expansion technique, we reduce the corresponding nonintegrable 
system of three coupled Gross-Pitaevskii equations (GPEs) to a completely integrable Yajima- 
Oikawa system. In this way, we obtain approximate solutions for small-amplitude vector solitons of 
dark-dark-bright and bright-bright-dark types, in terms of the mp = -f-l, —1, spinor components, 
respectively. By means of numerical simulations of the full GPE system, we demonstrate that these 
states indeed feature soliton properties, i.e., they propagate undistorted and undergo quasi-elastic 
collisions. It is also shown that, in the presence of a parabolic trap of strength uj, the bright 
component(s) is (are) guided by the dark one(s), and, as a result, the small-amplitude vector soliton 
as a whole performs harmonic oscillations of frequency ujjypl in the shallow soliton limit. We 
investigate numerically deviations from this prediction, as the depth of the solitons is increased, as 
well as when the strength of the spin-dependent interaction is modified. 



I. INTRODUCTION 

The development of far-off-resonant optical techniques 
for trapping of ultracold atomic gases has opened new 
directions in the studies of Bose-Einstein condensates 
(BECs) , allowing one to confine atoms regardless of their 
spin (hyperfine) state, see, e.g., Ref. One of major 
achievements in this direction was the experimental cre- 
ation of spinor BECs [2, 3], in which the spin degree of 
freedom (frozen in magnetic traps) comes into play. This 
gave rise to the observation of various phenomena that 
are not present in single-component BECs, including for- 
mation of spin domains 3 and spin textures [5|. 

A spinor condensate formed by atoms with spin F is 
described by a {2F + l)-component macroscopic wave 
function. Accordingly, a number of theoretical works 
have been dealing with multi-component (vector) soli- 
tons in = 1 spinor BECs. Bright [1,0, H and dark 
solitons, as well as gap solitons [10|, have been predicted 
in this context (the latter type requires the presence of 
an optical lattice). However, mixed vector soliton solu- 
tions (in particular, ones composed of bright and dark 
components) for the respective system of coupled Gross- 
Pitaevskii equations (GPEs) have not been reported yet, 
to the best of our knowledge. Actually, compound soli- 
tons of the mixed type may be of particular interest, as 
they would provide for the possibility of all-matter-wave 
waveguiding, with the dark soliton component forming an 
effective guide for the bright component, similar to the 
all-optical waveguiding studied in detail (chiefly, theo- 
retically) in nonlinear optics [llj . Such waveguides could 
be useful for applications, such as quantum switches and 
splitters emulating their optical counterparts p^ . 

On the other hand, mixed solitons of the dark-bright 
type were considered in a model of a two-component con- 
densate, described by two coupled GPEs {l^. Actually, 



the model also assumed that the two components repre- 
sented different spin states of the same atomic species, 
with equal scattering lengths of the intra-component and 
inter-component atomic collisions (i.e., the matrix of the 
nonlinear coefficients in the coupled GPEs was of the 
Manakov's type, that makes the system integrable in the 
absence of the external potentials). 

In this work we consider a quasi-one-dimensional 
(quasi-lD) F — 1 spinor condensate, described by a sys- 
tem of three coupled GPEs. In the physically relevant 
case of ^^Rb and ^"^Na atoms with F = 1, which are 
known to form spinor condensates of ferromagnetic and 
polar types, respectively (the definitions are given be- 
low), the system includes a naturally occurring small pa- 
rameter, S, namely, the ratio of the strengths of the spin- 
dependent and spin-independent interatomic interactions 
[ij, [l^. In the case of S — and without the exter- 
nal potential, the system of the three coupled GPEs be- 
comes the so-called Manakov system fl^, which is com- 
pletely integrable Exploiting the smallness of S, we 
will develop a multiscale-expansion method to asymp- 
totically reduce the nonintegrable GPE system to an- 
other completely integrable one, viz., the Yajima-Oikawa 
(YO) system. The latter was originally derived to de- 
scribe the interaction of Langmuir and sound waves in 
plasmas [lH and has been used in studies of vector soli- 
tons in the context of optics [ll| and binary BECs [20| . 
The asymptotic reduction is valid for homogeneous polar 
spinor BECs (such as ^'^Na), that are not subject to the 
modulational instability 0, IMI- Borrowing exact soli- 
ton solutions of the YO system, we predict two types of 
vector-soliton complexes in the spinor condensate, viz., 
dark-dark-bright (DDB) and bright-bright-dark (BBD) 
ones for the mp — -1-1,— 1,0 spin components. Numeri- 
cal simulations of the underlying (full) GPE system show 
that these solitary pulses (including ones with moder- 
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ate, rather than small amplitudes) emulate genuine soli- 
tons (solitary waves in integrable systems) quite well, 
propagating undistorted for long times and undergoing 
quasi-elastic collisions (quasi-elastic collisions of solitons 
in the nonintegrable two-component system were men- 
tioned earlier in Ref. [l3l|). 

The effect of the harmonic trapping potential on the 
solitons is also studied, analytically and numerically (the 
potential also breaks the exact integrability of the cou- 
pled GPE equations, even with (5 = 0). It is shown that 
the small-amplitude vector solitons of the mixed types 
perform harmonic oscillations in the presence of the trap 
of strength cj, at frequency a;/-\/2, which coincides with 
the well-known frequency of oscillations of a dark soli- 
ton in the single-component BEG In particular, the 
bright component (s) of the vector soliton oscillates at 
the same frequency, following their dark counterpart (s) 
(ordinary bright solitons in single-component BECs os- 
cillate at frequency ui, rather than |23j], accord- 
ing to the Kohn's theorem [Hj). As a matter of fact, 
it is a manifestation of the guidance of the bright com- 
ponent by the dark one. The oscillations of moderate- 
and large-amplitude solitons are studied as well and we 
find that the oscillation frequency is down-shifted from 
its value ui/V^ as the soliton depth is increased. In the 
particular case of large-amplitude solitons (that perform 
small- amplitude oscillations around the trap center), the 
frequency down-shift can be very well approximated by 
the prediction of Ref. [l^ . 

We also investigate the effect of a larger normalized 
spin-dependent interaction strength 6 on the stability of 
the predicted vector soliton solutions. In fact, we con- 
sider a value of 5 an order of magnitude larger than the 
physically relevant one of the polar sodium spinor con- 
densate {S = 0.2 rather than 0.0314) and investigate nu- 
merically the respective coupled GPEs. The result is 
that, generally, under such a strong perturbation the 
solitons emit stronger radiation and are eventually de- 
stroyed. However, even in such a case of a large 6, small- 
and moderate-amplitude DDB solitons are found to per- 
sist up to relatively large times, of order of more than 
300 ms in physical units. Given that for the physically 
relevant small value of S pertaining to the sodium spinor 
BEG, the respective lifetime was four times as large, we 
believe that the vector solitons derived in this work have 
a good chance to be observed experimentally. 

The paper is organized as follows: In Sec.[II]we present 
the model, expound our analytical approach for the ho- 
mogeneous system, and derive solutions for the bright- 
dark soliton complexes. Section|ni]is devoted to the pre- 
sentation of the numerical and analytical results for the 
dynamics of the solitons in both the homogeneous and 
inhomogeneous (harmonically confined) system. Finally, 
Sec. IIVI concludes the paper. 



II. THE MODEL AND ITS ANALYTICAL 
CONSIDERATION 

A. The model 

At sufficiently low temperatures (finite temperature ef- 
fects have been considered recently in Refs. [1^), and 
in the framework of the mean-field approach, the spinor 
BEG with = 1 is described by the vector order param- 
eter, *(r,t) = [1'_i(r, i),*o(i",i),*+i(i", i)]^, with the 
components corresponding to the three values of the ver- 
tical spin projection, mp = —1, 0, -1-1. Assuming that the 
condensate is confined in a highly anisotropic trap with 
frequencies lo^ ^ , we may assume the wave functions 
approximately separable, ^'o,±i ~ ''pQ,±iix)'ip±{y, z), 
where the transverse wave function ip±iy, z) is the ground 
state of the respective harmonic oscillator. Then, av- 
eraging of the underlying system of the coupled three- 
dimensional (3D) GPEs in transverse plane (y, z) [2^ 
leads to the following system of coupled ID equations for 
the longitudinal components of the wave functions (see 
also Refs. 0, 0, B i, IS] ) : 

+ crVo'^;i, (1) 

+ 24i^V-iV'o^+i, (2) 

where the asterisk denotes the complex conjugate, and 
H^i = -{ti^/2m)dl + {\/2)mulx^ + c[,^^^ntot is the 
spin-independent part of the Hamiltonian, with ritot = 
IV'-iP + IV'oP + IV' -Hip being the total density (m is 
the atomic mass). The nonlinearity coefficients have 

an effectively ID form, namely c'q^^ = co/27ra^ and 

c^2^^ = C2/27ra^, where = yJJ^JrmJZ is the transverse 
harmonic oscillator length, which defines the size of the 
transverse ground state. Finally, coupling constants Cq 
and C2 which account, respectively, for spin-independent 
and spin-dependent collisions between identical spin-1 
bosons, are given by (in the mean-field approximation) 

Co = 5 , C2 = , (3) 

3m 3m 

where oo and 02 are the s-wave scattering lengths in the 
symmetric channels with total spin of the colliding atoms 
F = and F = 2, respectively. Note that the F 1 
spinor condensate may be either ferromagnetic (such as 
the *^Rb), characterized by C2 < 0, or polar (such as the 
23Na), with C2 > [13,123. 

Measuring time, length and density in units of 

h/c^Q^'^no, h/\J mc^^^'no and no, respectively (where no 
is the peak density), we cast Eqs. in the dimen- 
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sionless form, 

+ si^^or^i, (4) 
idt4>o = FsiV'o+^dV'-ii'+iv^+inv'o 

+ 25^_iV'o*^+i, (5) 

where i?si = — (l/2)c^x+(l/2)f^tr2^^+"-tot, the normahzed 
harmonic trap strength is given by 



also be stressed that, although stationary equations pre- 
sented below may indeed be found from the system of 
two, rather than three, coupled GPEs, the stability tests 
for the solutions are performed against general perturba- 
tions (see below), which include those with -0+1 ^ ip-i 
(i.e., the full system of the three equations was employed 
in the stability simulations). 



B. Linear analysis 



w_L 2(ao -I- 2a2)no ' 



and we define 



0-2 - QQ 

ao -I- 2a2 



(6) 



(7) 



According to what was said above, 5 <Q and (5 > cor- 
respond, respectively, to ferromagnetic and polar spinor 
BECs. In the relevant cases of ®^Rb and ^"^Na atoms with 
F =1,6 = -4.66 X 10-3 [13 and (5 = -f3.14 x 10"^ 
respectively, i.e., in either case, 5 plays the role of a small 
parameter of Eqs. (l4])-([5]). 

Generally, Eqs. (H))-® give rise to spin-mixing states 
[29'! . However, there also exist non-spin-mixing, or spin- 
polarizcd states, which are stable stationary solutions of 
Eqs. (HI) and © [13, HO]. Here, we first consider the 
spatially homogeneous system (fitr = 0), and focus on 
such solutions having at least one component equal to 
zero, the remaining ones being continuous waves (CWs). 
The corresponding exact stationary solutions are 



and 



ip^i = V+i = Y 2 '^^P(^*Ai^), fPo = 0, (8) 
■0-1 =-0+1=0, 00 = VAiexp(-i/it). (9) 



As we demonstrate below, weakly nonlinear perturba- 
tions around these solutions take the form of three- 
component dark-bright soliton complexes. In particular, 
perturbations around solutions ^ and ([9|) will lead to 
solitons of the DDE and BED types, respectively, for 
components ■0±i and "00- Since the analytical approach 
and the derivation of the soliton solutions for the two 
cases are quite similar, we focus below on the DDB soli- 
tons, and discuss the BED ones in a brief form. 

It is relevant to note that, for solutions with 0+i = 
tA-i = -01 J two equations ^ coalesce into a single 
one for wave function "01. Then, the transformation 

V'l = (0D + 0b) / {2VTTP), Vo = (0D - 0b) /VT+P 

casts the system of two equations ^ and ([5]) into the 
form of two coupled GPEs, with nonlinear cross-coupling 
coefficients go =93 = (l~<5)/(l + '5)i which was intro- 
duced in Ref. [l^ (also with the objective to study bound 
complexes of dark-bright solitons, carried by fields 0n 
and 0B, respectively). In fact, the analysis in Ref. [l3[ 
was limited to the case oi go = Sb = 1, while in this work 
we focus on effects generated by small J 7^ 0. It should 



Aiming to find solutions of Eqs. (HI)-® close to the 
CW solution given by Eq. ([5]), we start the analysis by 
adopting the following ansatz. 



^+1 
■00 



\/ n{x, t) exp[— i/it -f- i0(2;, t)], 
<i>o(a;, i) exp(-i/it), (10) 



where n{x,t) and 4>{x,t) are real density and phase of 
fields '0±i, while $0 is, generally, a complex function. 
Substituting Eq. ^ into Eqs. O-®, we derive the 
following a system. 



;[dtn + d^ind^cj))] - n [9^0 + 2n - /i + (1 + 5)\^af 



= 0, 



-a2$o-(27i-Ai+|<i>on<i>o 
26n ($0 + ^e"^"!') = 0. 



(11) 



(12) 



The CW state ([8]) corresponds to an obvious solution of 
Eqs. (HU and ^ with n = /x/2, = 0, $0 = 0. Next, 
we linearize the equations around this state, looking for 
a solution as n = (/^/2) + en,, = e0 and $0 = e$o, 
where e is a formal small parameter. At order 0(e), the 
linearization leads to the following system: 



i (dtn + |a20) - (dt^ + 2fi- ^dln^ = 0, (13) 



$S =0. 



(14) 



Combining real and imaginary parts of Eq. (|13p . we arrive 
at a dispersive wave equation. 



(15) 



which gives rise to a stable dispersion relation between 
wavenumber k and frequency lu (the absence of complex 
roots for uj at real k implies the modulational stability of 
the underlying CW state): 



(16) 



It follows from Eq. that, in the long- wave limit (fc 
0), small-amplitude waves can propagate on top of CW 
solution ([8|) with the speed of sound, 



(17) 
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A similar analysis for Eq. (fTi]) . which is decoupled from 
Eq. (fT5|) . leads to dispersion relation 



(18) 



It is clear from here that, for (5 > (which corresponds 
to the polar state), Eq. (fT8)) has no complex roots for w, 
hence the trivial solution to Eq. (fT2| . (f>o = 0, is modu- 
lationally stable. However, i5 < (corresponding to the 
ferromagnetic state) gives rise to modulational instability 
of the $0 = solution against the perturbations whose 
wavenumbers belong to the instability band, k < 2yJ\5\^. 
Note that these results comply with those reported in 
Ref. [llj. Below, we focus on the modulationally stable 
case, which pertains to the polar state with 5 > Q. 



C. Asymptotic soliton solutions 

To consider solutions for weakly nonlinear deviations 
from the CW state, we recall that 5 is a small param- 
eter of the system (|1])-(IS]), which suggests to define the 
stretched variables, 

X = 5^/^{x-^it), T=6t. (19) 

Then, we seek for solutions of Eqs. (fTT|l -(fT2 |) as 

n = (/i/2) + J • p, (f> = S^^^a, ^Q^S^'^^q, (20) 



q = qi cos{Kx ~ J7i) + iqi sin(_R'x — fJi), 



(21) 



where p = p(A, T), a = a(A, T), gi_2 = q\^i(X, T), while 
K and fl are unknown wavenumber and frequency. Sub- 
stituting Eqs. ((20|) in Eq. (fTTj) . at the leading order in 6, 
which is 0((5), we derive a relation between density p and 
phase a, 



y/JIdxa = 2p. 



(22) 



At the next order 0{S'^^^), the resulting equation is com- 
plex: 

- iip/4) {2dxp ~ ^d\a) + dra + \q\^ - 0. (23) 

The imaginary part of the expression on the left-side of 
Eq. ((23)) vanishes due to the validity of Eq. ((22)) . while 
the real part leads to equation 9tQ! + \q^ = 0. The 
compatibility condition of the latter with Eq. (I22p leads 
to 



dTP=-{^l2)dx{\q\ 



(24) 



We now proceed to Eq. ((T^), which, to leading order 
in (5, i.e., at 0{8^l'^), yields the following system: 



riqi - (if 2/2) 92 = 0, 



Nontrivial solutions to Eqs. (|25p are possible when the 
following dispersion relation for and K holds: 

V? ^ {p + K^ ji) . (26) 

Next, to order 0{6^l^), Eq. HH) leads to system 

- + Kdxq2 = 0, 

-Kdxqi + ^/pdxqi = 0, 

which has nontrivial solutions if = p. In combination 
with Eq. (|26p . the latter relation selects the frequency, 
n = 5^74- Finally, at order 0(5^/^), Eq. ^ leads to 
equation 



(27) 



Equations ([24]) and ((27)) . which are the basic result 
of our analysis, constitute the Yajima-Oikawa (YO) sys- 
tem. It describes the interaction of low-frequency and 
high-frequency waves, and was originally derived in the 
context of plasma physics, where it applies to Lang- 
muir (high-frequency) waves, forming a packet (soliton) 
moving at velocities close to the speed of sound, and 
thus strongly coupled to the ion-acoustic (low-frequency) 
waves [ll]. As shown in Ref. [ll|, the YO system is 
integrable by means of the inverse-scattering-transform, 
and gives rise to soliton solutions. The solitons have the 
— sech^ shape for field p, and sech shape for q, which cor- 
respond to a density dip for components ijj±i and a bright 
soliton for ipo, as per Eqs. ([20]) . According to Eq. (|22p. 
the phase profile of the '>p±i components, in the form of 
tanh, is associated to the density dip, hence the patterns 
in these components, generated by the exact solution of 
the YO system, are genuine dark solitons. The full form 
of the approximate (asymptotic) DDB soliton solution 
to Eqs. (HI)-®, into which the YO soliton is mapped by 
Eqs. (Uni), (Uni), and (EOl), is 



^±iix,t) = \/{p/2) - 2STj'^sech^{2V5TjZ) 

X exp —ipt— 2iT]^/ 5 / pta.nh{2V6r]Z) , (28) 

Mx,t) = 2=^/2^=^/4^^-i/4y^gg^j^^2V^r;Z) 



X exp 



-ipt + i^x - 2iVd^Z + 2i5{rj^ - ^'^)t \ , (29) 



(25) 



where Z = x — (y/p — 2y/5S^)t, while 77 and ^ are arbitrary 
parameters of order 0(1). 

A similar analysis can be performed to derive asymp- 
totic soliton solutions of the BED type. In that case, 
starting from CW solution ([H]), we seek for solutions of 
Eqs. ©-(in]) in the form of 

V'-i = V'+i = *o(a;,t)exp(-ipt), 

= \/ n{x, t) exp[—ipt + i(l){x, t)]. (30) 
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Next, following the same analytical approach which has 
led above to the DDB soliton, we again end up with 
the YO system, in a form similar to Eqs. ([24|) and ([27|) . 
namely, 



idrq + \d\q - pq = 0. 



(31) 



Eventually, the approximate BBD soliton solution to 
Eqs. (HI)-®, generated by the YO soliton, is 



X exp 



~ipt + iy/JIx - 2iVS^Z + 2i5{Tf ~ C^)tj , (32) 



ipo{x,t) = \J {p/2) - 4(57/2sech2(2%/^77Z) 
-ipt — 2ir\\/ 5 / p tanh(2'\/5r/Z) 



X exp 



(33) 



As the latter solution is quite similar to the DDB one, 
given by Eqs. (p8)) - (|29|) . below we only deal with the dy- 
namics of the DDB solitons. It is worthwhile to note in 
passing that both types of solutions are genuinely travel- 
ing ones, i.e., they do not exist with zero speed. 



III. DYNAMICS OF DDB SPINOR SOLITONS 

A. Numerical Results 

In order to test the prediction of the existence of the 
DDB solitons in the underlying spinor BEG model, we 
turn to numerical integration of the original GPEs 
([5]). In particular, we fix 5 = 0.0314 (the value corre- 
sponding to ^■^Na) and use the following initial conditions 
for the densities, 

1 



|V'±i(x,i = 0)p = -[p-u&ec\-?{^x)\, (34) 
\il)o{x,t = Q)\^ = -prsecYi'iy/Dx), 



(35) 



while initial phase profiles are similar to those in 
Eqs. the parameter that determines the ini- 

tial width of the soliton being 



4r/2,5. 



(36) 



Other parameters are chosen as = 2, ^ = 0.5, and 
f^tr = (for the homogeneous condensate), or fitr = 0.05 
for the trapped condensate. In physical terms, this choice 
corresponds to the spinor condensate of sodium atoms 
with the peak ID density no ~ 10^ m^^, which con- 
tains ~ 20000 atoms confined in the trap with frequen- 
cies L0± = "i^LOx = 27r X 230 Hz; in this case, the time and 
space units are, respectively, 1.2 ms and 1.8 /xm. 



Choosing value 77 = 1 of the arbitrary parameter intro- 
duced above, and substituting 5 = 0.0314, u = 0.13, we 
have checked that these values indeed provide for very 
good agreement of the analytical predictions with nu- 
merical results. However, we will display numerical re- 
sults obtained for an essentially larger value of viz., 
V = 1.2 [which, as seen from Eq. ([M]) corresponds to 
77 = 3.091 and hence, from Eq. to a soliton com- 

plex with deeper and narrower dark components, and 
similarly taller and narrower bright components] . In this 
way, we intend to showcase the really wide range of va- 
lidity of the analytical approach, and the robustness of 
the obtained solitary wave solutions. 

More specifically, we first check if the spinor DDB soli- 
ton complexes indeed behave as genuine solitons, in the 
small-amplitude limit. To this end, we take initial con- 
ditions in the form of the superposition of two different 
pulses, 



X exp { —i\ — tanh(a::+) + i, — tanh(a;_) ) , (37) 
\ M p V M / 



sech(x+)e^ 



-hsech(a;_)e-'V¥^-+'(«/'')^- 



(38) 



where x± = v^(x±a;o) and xq = ±15 are positions 
of centers of the two pulses. As seen in Eqs. ([57)1 - ([55)1 . 
the soliton components are lent opposite initial momenta 
and, as a result, they propagate in opposite directions, 
as shown in the top panel of Fig. [TJ We stress that 
even though a small amount of radiation is emitted (see 
the four bottom panels of Fig. [T|) , the two dark solitons 
in the V'ii fields, coupled to their bright counterparts 
in the V-'o component, propagate practically undistorted, 
and around t — 19 they undergo a quasi-elastic collision; 
moreover, it is clearly observed that the solitons remain 
unscathed after the collision. This result is consistent 
with our asymptotic calculations performed above, in- 
dicating that the small-amplitude limit (for sufficiently 
small S) of the nonintegrable system of Eqs. (HI)-® be- 
haves like the integrable YO system. 

Next, we consider the confined system, with f2tr = 
0.05. In this case, strictly speaking, the asymptotic re- 
duction of the primary system of Eqs. H])-® to the YO 
system [Eqs. ((M)) and P7)) ] is not valid. However, even in 
the presence of the potential term, the solutions obtained 
with r^tr = may be used as an initial configuration set 
near the bottom of the trap, to generate DDB- (or BBD-) 
like solutions of the inhomogeneous system. To that end, 
we first integrate Eqs. ®-® in imaginary time, finding 
a ground state of the Thomas-Fermi (TF) type for fields 
'0±i, which is approximated by the well-known analyti- 
cal density profile [31], n±i = {l/2){fi — n^j.x^). Then, at 
t = 0, the initial conditions for the ijj±i components are 
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FIG. 1: The two top panels show contour plots of the densi- 
ties of the tp±i (left panel) and ^po (right panel) components of 
the spinor condensate (with 5 = 0.0314) in the homogeneous 
system (f2tr = 0). The tp±i components contain a pair of dark 
pulses, initially placed at xo — ±15, that, together with the 
bright components in the ipo field coupled to them, undergo 
a quasi-elastic collision at t ~ 19 and propagate unscathed 
afterward. The parameters are jj, = 2, ^ = 1.54, rj = 3.091, 
and V — 1.2. The four bottom panels show snapshots of the 
densities observed at t — 0, 10 (before the collision), t = 19 
(when the collision occurs) and 40 (after the collision). 




FIG. 2: The two top panels show contour plots of the den- 
sities of the tf}±i (left) and t/jo (right) fields confined in the 
harmonic trap with Qti = 0.05 (the other parameters are the 
same as in Fig. [l|. Initially, each of the Thomas- Fermi pro- 
files of the tp±i components carries a dark soliton, while the 
i/iQ component is a bright soliton (the initial position is at the 
trap's center, a; = 0). The four bottom panels show snapshots 
of the densities observed at t = 0, 444, 937, and 987. 



B. The soliton's oscillation frequency 



taken as the numerically found TF profiles multiplied by 
the dark soliton as in Eq. while the initial config- 

uration of the V'o field is taken as the bright soliton in 
Eq. 

In such a case, and given that the spinor DDB solitons 
were found above to be robust objects behaving similar 
to genuine solitons of an integrable system, one would 
expect that the solitons would perform harmonic oscilla- 
tions in the presence of the (sufficiently weak) parabolic 
trap. Figure [2] shows that this is the case indeed: The 
DDB soliton, which was initially placed at the trap's cen- 
ter, oscillates as a whole without significant deformations 
of its components up to large times [while the figure ex- 
tends to t = 1000 (which is 1.2 seconds in physical units), 
a similar behavior continues at still larger times]. This 
is a clear indication to the fact that the predicted DDB 
complexes have a good chance to be observed in the ex- 
periment. A noteworthy feature of the numerical data is 
that the bright-soliton component is guided by the dark 
ones, the entire soliton complex oscillating at a single fre- 
quency. The value of the frequency is estimated below. 



The effect of the harmonic trap on the spinor-soliton 
dynamics can be studied analytically, using the asymp- 
totic multiscale expansion, similar to how it was done in 
Refs. [s^jIH, H^l- In those works, asymptotic reductions 
of various scalar GPE-based models, which included the 
trapping potential, led to Korteweg-de Vries equations 
with variable coefficients [instead of constant ones that 
would be the case for the homogeneous (untrapped) sys- 
tem]. In the present situation, we may expect, accord- 
ingly, that the inclusion of the harmonic potential may 
lead to a YO system with variable coefficients. 

An important result that can be produced by such an 
analysis is the oscillation frequency of the solitons in the 
presence of the harmonic trap. As shown in Refs. (3^ . 
i33i . i34j b y m eans of the adiabatic perturbation theory for 
solitons [23], the oscillation frequency can be obtained 
from the equation (which is valid to leading order in the 
small parameter characterizing the inhomogeneity of the 
system), 

§-HX). (39) 

where X is a properly chosen slow spatial variable, and 
c is the speed of sound, which is now position-dependent 
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due to the background density of the trapped condensate. 
In fact, Eq. (|39p is a straightforward generahzation of 
the result obtained for the homogeneous system, where 
it has been found that the speed of sound is c = ^/I [see 
Eq. P?)) ]. while the soliton's velocity is u = ^/Ji — 2%/^^, 
implying that (for sufficiently small S) dxj dt ~ ^JJi = c. 

To adopt this approach to the present case, we need 
to find the local speed of sound c{X) when the har- 
monic potential term, V = {l/2)U'^^x^ , is included in 
the spin- independent part of the Hamiltonian, Hsi- Tak- 
ing into regard that the potential varies slowly on the 
soliton's spatial scale, (see, e.g.. Fig. [5]), we define 

the above-mentioned slow spatial variable as X = ex, 
where e = iltr/i^tr [recall that ilti given in Eq. ([6|) is of 
order 10~^], and fttr is an auxiliary 0(1) scale parame- 
ter. This way, the trapping potential takes the form of 
V{X) = (l/2)J7^j.X^, i.e., it depends only on slow vari- 
able X. 

Then, the local speed of sound can easily be derived 
upon considering the linearization of Eqs. pT|) and ([H]), 
which are modified by the inclusion of the term —nV{X) 
in the left-hand side of Eq. pT|) . The ground state of 
this system can easily be found by setting v = 4'x = 
and (j^t = —fJ-- Then, since Eq. (|lip implies that n — hq 
is time-independent in the ground state, we assume that 
no — no{X) and, to the leading order in e, we obtain 



no(^) = (l/2) f.i~ViX) 



(40) 



in the region where /j > V{X), and uq — outside. 
Equation ([3D]), which is the TF approximation for the 
density profile, also impHes that, for V{X) = (l/2)fij,.X^ 
the axial size of the trapped condensate is 2L = 2-^2/7/0. 
Similarly to the analysis presented above in Sec. IIIBl 
we now consider the linearization around the ground 
state and seek for respective solutions to Eqs. ([TT|) - (fT^ 
as n = no(X) + ini{x,t), (j) = — A*oi + e<?!'i(a^, t), and 
$0 = e$i(a;,t), with ~ exp[i(fcx - wt)]. This 

way, we obtain the following dispersion relation for the 
inhomogeneous system. 



^2 ^ 2no{X)k^ + (l/4)fc'', 
and, accordingly, the local speed of sound: 



z{X) = ^j2no{X), 



(41) 



(42) 



which bears resemblance to the sound propagation in 
weakly nonuniform media ^36] ; in the homogeneous case, 
Eq. ^ is reduced to Eq. (fTT]) . 

Next, we substitute Eq. (gl]) in Eq. and, taking 
into regard the density profile given by Eq. ((40l) , we inte- 
grate the resulting first-order differential equation. The 
result is 




0.5 
V/|J. 



FIG. 3; Left panel: Contour plot of the density of a shal- 
low dark soliton (shown are the identical ji/Jiij^ fields) with 
initial amplitude v = 0.13; the other parameter values are 
77 = 1, ^ = 0.5, while the values of the chemical potential 
and the harmonic trap strength are jj, — 2 and Qtr = 0.05 
as before. The soliton performs oscillations with a frequency 
i^osc ~ 0.03513, which is almost identical to the analytical 
prediction of ntr/V2 = 0.03536 (error ^ 0.65%). Right 
panel: Percentage error of the numerically found soliton os- 
cillation frequency, as compared to the analytical prediction 
of Eq. (|44|l . as a function of the relative dark soliton depth 
v/l-t- The region ly/fi ^ 1 corresponds to shallow solitons, 
while v/n = 1 corresponds to a black soliton; for the case 
shown in the left panel u/fi — 0.065. As seen, in the for- 
mer case the analytical prediction is fairly good (the error is 
below 5% for every v/fi < 0.25), while for large amplitude 
solitons p/fi > 0.7 the error becomes more than 10%, with a 
maximum 13.7% for v/jj. = 1. 



which demonstrates that a sufRciently shallow spinor 
dark soliton oscillates with frequency 



71' 



(44) 



X = L sin 



(f]tr/V2) t 



(43) 



similarly to the known result for the oscillations of dark 
solitons in a single-component BEC [l^l [note that, for 
the sake of simplicity, we dropped the tilde in Eq. (gH)]. 
We should clearly stress here that we expect the above 
result to be valid in the case of shallow solitons, with 
a relative depth i//^ <^ 1 [see Eq. ([M]) ]. For example, 
as seen in the left panel of Fig. [31 in the case of such a 
shallow soliton with v/ji ^ 0.065 (corresponding to the 
above mentioned physically relevant choice of = 0.13 
and a value of the chemical potential fi = 2) the analyt- 
ical prediction is quite accurate: the numerically found 
oscillation frequency (for fitr = 0.05) is « 0.03513, while 
the analytical prediction of Eq. is 0.03536 (the per- 
centage error is « 0.65%). On the other hand, as seen 
in the same figure, the amplitude of the soliton oscilla- 
tion is 39.992, while L = y^%l/n = 40 (the respective 
percentage error is 0.02%. 

Although the result of Eq. ([33]) [and ([33])] is "com- 
patible" with our analytical approach developed in the 
previous section, which is also valid for small-amplitude 
solitons, it is also of interest to numerically investigate 
the oscillations of solitons of moderate and large ampli- 
tudes. In this respect, in the right panel of Fig.[3]we show 
the percentage error of the numerically found soliton os- 
cillation frequency [relative to the prediction of Eq. (|44p ] 
as a function of the relative dark soliton depth v/fi. It 
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is clear that the region v/ <Si \ corresponds to shallow 
solitons, while the limiting case oi v / ^ = 1 corresponds 
to a "black" soliton, with an initially zero intensity at 
the trap's center (the latter is slightly displaced from the 
trap's center to be set into motion). As seen in the fig- 
ure, the analytical prediction of Eq. (I44|) is very good as 
concerns the estimation of the oscillation frequency, since 
the error is below 2% for every v/ ^ < 0.1. Notice that we 
have also computed the error in the oscillation amplitude 
(not shown here), which we have found to be larger (i.e., 
up to 8% in this regime); this is due to the fact that the 
increasingly deeper solitons are not reflected at the rims 
of the condensate, but rather inside the cloud, as can be 
seen e.g., in Fig. [2l 

In the case of moderate (and large) amplitude soli- 
tons the analytical prediction becomes worse: For exam- 
ple, in the particular case shown in Fig. [2] (in this case 
v/fi — 0.6), comparing the numerically found soliton's 
oscillation frequency, Wosc ~ 0.032 to the aforementioned 
prediction of 0.03536 (again for Qtr = 0.05), we find a rel- 
atively large discrepancy (« 9.5%) between the two val- 
ues. However, an important observation regarding Fig. [21 
is that the bright-soliton component performs oscillations 
at the same frequency, ^tr/V^. This is a clear indication 
that the bright component is guided (effectively trapped) 
by the dark component of the DDB complex. Note that, 
in the single-component BEC, bright solitons oscillate 
in the parabolic potential with frequency, iltr [23j (as 
a consequence of the Kohn's theorem [2J|), rather than 
nu/V2. 

Furthermore, we observe that, naturally, the discrep- 
ancy becomes even larger in the case of large-amplitude 
(nearly-black) solitons, which perform small-amplitude 
oscillations around the trap's center. For example, for 
v/fi = 0.8 {l'/^J- = 0.9) the numerically found values 
of the oscillation frequency deviate from the value of 
ntr/V2 by 11.5% (12.4%), while in the limiting case of 
v/fx = 1 the respective percentage error is 13.7%. It is 
clear that such large deviations are due to the fact that 
the above numerical results pertain to non-small ampli- 
tude solitons with large values of fJ.^ contrary to what 
is the case for the the analytical approach which requires 
u/fi <ti 1, as mentioned above. 

In such a case of large amplitude solitons, it is inter- 
esting to compare the numerically found oscillation fre- 
quencies to a different analytical prediction presented in 
Ref. [l^ . In this work, the oscillation frequency may re- 
sult from a nonlinear equation of motion for the bright- 
dark soliton in a binary BEC mixture (see Eq. (5) of 
[l3|). In fact, for shallow solitons the oscillation fre- 
quency is identical to the one given by Eq. (|44p . while, 
in the opposite limit of very deep dark solitons, may be 
approximated (in our language) as, 

where A'b is the number of atoms of the bright-soliton 



component, which, employing Eq. (|35p . is easily found to 
be A^B = 2i/^/^^/?7^. It is clear that Eq. (j45]) shows that 
the oscillation frequency is down-shifted (as compared to 
the value of flti/V^) [i-e., the dark-bright pair executes 
slower oscillations, as the bright component is enhanced], 
thus having a quantitative agreement with our numerical 
observations. Perhaps more importantly, it also provides 
better quantitative estimates for the values of the soli- 
ton's oscillation frequencies. In particular, for normalized 
soliton depths v/fj, = 0.8, 0.9 and 1, the predictions of 
Eq. (|45p deviate from the respective numerically found 
oscillation frequencies with percentage errors 1.7%, 4% 
and 5.2% (recall that a respective comparison with the 
prediction of Eq. (|^ led to percentage errors 11.5%, 
12.4%, and 13.7%, respectively). The above results in- 
dicate that a more detailed description of the motion of 
the bright-dark soliton complexes in the trapped spinor 
condensate, in the lines of Ref. [l3j , would be very inter- 
esting. However, this is beyond the scope of the present 
work. 



C. Effects of larger spin-dependent interactions 

In the previous subsections we had kept the value of 
the parameter S fixed, i.e., S = 0.0314 (corresponding 
to the polar spin-1 sodium BEC). Such a small value of 
6 validates our perturbative approach, which allowed us 
to find approximate DDB soliton solutions of the YO- 
type, and study their oscillations in the trapped spinor 
condensate. It is interesting, however, to test the stability 
and the dynamics of the obtained DDB solitons in the 
more general case of a stronger perturbation, considering 
also non-small values of S. In this respect, we will now 
present numerical results obtained by direct numerical 
integration of Eqs. (j4]) and ([5]) for S = 0.2, which is an 
order of magnitude greater than the previous value. Note 
that we will study the evolution of DDB solitons with the 
same amplitudes as in the case of small S, so as to directly 
compare the results pertaining to weak or moderate spin- 
dependent interaction strength. 

In Fig. m we show the evolution of a moderate am- 
plitude DDB soliton, characterized by the parameters 
f — 0.61, rj = 1.22, and = 1.2; for the same val- 
ues of the trap strength and chemical potential as be- 
fore (Jltr = 0.05 and /i = 2), the initial condition is the 
same as the one shown in the second-row left panel of 
Fig. [2] (in this case, recall that the normalized amplitude 
of the dark solitons in the components is fi = 0.6). 
As observed in Fig. [H although the stronger perturba- 
tion induces emission of stronger radiation (see bottom 
left panel), the losses are not significant up to relatively 
large times, of order of t = 311 (or t — 360 ms in physi- 
cal units): the dark (bright) soliton densities are only 8% 
(9%) smaller than their initial values, and it can safely 
be concluded that even for such a strong value of 6 the 
DDB vector soliton has a good chance to be observed 
in an experiment (provided, of course, that such a mag- 
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FIG. 4: Same as Fig. [1 but for 5 = 0.2. The soliton pa- 
rameter values are ^ = 0.61, — 1.22, and i/ — 1.2, while 
the values of the trap strength and chemical potential are 
Jltr =0.05 and fj, = 2. This choice renders the initial soliton 
densities identical to the ones shown in the second-row left 
panel of Fig. [2] The two bottom panels show snapshots of 
the densities at t = 311 and t = 967. 



nitude of the spin-dependent interaction is experimen- 
tally achieved). However, for later times the continuous 
perturbation-induced emission of radiation results in the 
eventual destruction of the DDE complex. In particu- 
lar, at f « 1000, the solitons almost decay: at t = 967 
(see bottom right panel of Fig. d]), the density of the 
dark (bright) soliton is 63% (28%) smaller than its ini- 
tial value. 

The cases of large- and small-amplitude solitons, with 
normalized amplitudes (of the dark soliton) being v/ ii = 
0.8 and v/ fj, = 0.065, respectively, was examined too (re- 
sults not shown here). It was found that, naturally, the 
large-amplitude DDE soliton starts to accelerate immedi- 
ately due to the strong radiation emitted and decays fast, 
being completely destroyed at t sa 300: the soliton densi- 
ties are smaller than 50% of their initial values. On the 
other hand, the small-amplitude DDE soliton was found 
to be slightly more robust, showing a behavior similar 
to the one of the moderate amplitude soliton shown in 
Fig. m but for smaller times: in fact, the vector soliton 
persists up to t w 600 (as per the above criterion, the 
densities are smaller than 50% of their initial values) and 
then decays. 

To conclude this subsection, it is clear that the small- 
and moderate-amplitude YO-type DDE vector soliton 
persist in the spin-1 condensate up to experimentally 
relevant times even for large spin-dependent interaction 
strength, namely an order of magnitude larger than the 
physically relevant value pertaining to the polar sodium 
spinor condensate. 



IV. CONCLUSIONS 

We have studied bright-dark soliton complexes in po- 
lar spinor Eose-Einstein condensates, both analytically 
and numerically. Our analytical approach is based on 
the small amplitude-asymptotic reduction of the nonin- 
tegrable vector (three-component) system of the coupled 
Gross-Pitaevskii equations to a completely integrable 
model, viz., the Yajima-Oikawa (YO) system. Eorrow- 
ing soliton solutions of the YO system and inverting 
the reduction, we have obtained an analytical approxi- 
mation for small-amplitude vector solitons of the dark- 
dark-bright and bright-bright-dark types, in terms of the 
mp = +1,-1,0 components, respectively. The analyti- 
cal predictions were confirmed by direct numerical sim- 
ulations. The so constructed approximate soliton states 
were found to propagate undistorted and undergo quasi- 
elastic collisions, featuring properties of genuine solitons. 

The effect of the harmonic trapping potential (which 
also contributes toward the nonintegrability of the un- 
derlying equations) on the solitons was also studied nu- 
merically and analytically. It was found that even vec- 
tor solitons of moderate (non-small) amplitudes maintain 
their identity in the presence of the parabolic trap, and 
perform harmonic oscillations, at large times (~ 10 sec- 
onds or even more, in physical units). We have found 
that the soliton's oscillation frequency takes (in the an- 
alytical approximation) the value of flti/V^, the same 
as featured by the dark soliton in the single-component 
repulsive trapped condensate. The numerical results ver- 
ify the analytical prediction, for sufficiently shallow dark 
solitons: It was found that for initial soliton depths be- 
low 10% of the chemical potential, the deviation of the 
analytical prediction from the numerically found oscilla- 
tion frequency was below 2% (the error in the estimation 
of the amplitude of the soliton oscillation was below 8%). 
For moderate and large amplitude solitons, the discrep- 
ancy in the frequency was larger (of the order of more 
than 10%); however, in the case of very deep dark soli- 
tons we found that the respective prediction of Ref. [l^ 
led to a significantly smaller error, less than 5%. Thus, 
an elaborated description of the bright-dark soliton mo- 
tion (for solitons of arbitrary amplitude) in the trapped 
spinor condensate is certainly a challenge for future work. 

We also tested the robustness of the derived vector 
soliton solutions in the case of a large normalized spin- 
dependent interaction strength, namely an order of mag- 
nitude larger than the physically relevant one correspond- 
ing to the polar sodium spinor condensate. We found 
that although the solitons are generically destroyed un- 
der such a strong perturbation, the lifetimes of small- 
and moderate-amplitude DDE solitons was more than 
300 ms in physical units. Thus, the vector solitons ob- 
tained in this work have a good chance to be observed in 
an experiment. 

The bright-soliton component (s) was found to be 
guided by their dark counterpart (s), oscillating with the 
frequency imposed by dark components. This is an ex- 
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ample of the all-matter-wave soliton guidance, with po- 
tential applications in the design of quantum switches 
and splitters. 
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